Phyllotactic description of hard sphere packing in cylindrical channels 
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We develop a simple analytical theory that relates dense sphere packings in a cylinder to cor- 
responding disk packings on its surface. It applies for ratios R = D/d (where d and D are the 
diameters of the hard spheres and the bounding cylinder, respectively) up to _R = 1 + l/sin(7r/5). 
Within this range the densest packings are such that all spheres are in contact with the cylindrical 
boundary. The detailed results elucidate extensive numerical simulations by ourselves and others 
by identifying the nature of various competing phases. 
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Packing problems, such as that of finding the dens- 
est arrangement of hard spheres, are among the most 
ancient and challenging problems in mathematics and 
science. The most famous example is the Kepler con- 
jecture (1611), concerned with identical spheres in an 
unbounded volume p]|. In contrast, the optimal packing 
of such spheres within a hounded space is much less well 
understood. A prominent example is finding the densest 
packing arrangement of monodisperse spheres within a 
cylindrical channel, which we consider here. There are 
extensive results from simulations, which we corroborate 
and augment in this Letter. We supply a detailed in- 
terpretation of them, by making a connection with the 
analogous (but simpler) problem of packing disks on the 
bounding surface, which is an interesting problem in its 
own right 

Our procedure is therefore to first explore an analytical 
treatment of the surface disk packing problem, and then 
develop a transformation which gives approximate but 
illuminating results for sphere packing within a cylinder. 
Such surface patterns have their origin in the early study 
of intriguing spiral arrangements of leaves and stems [Sj . 
Thus they are associated with the term phyllotaxis which 
now stands for the standard classification of triangular 
patterns on a cylinder, in any context. It is defined and 
used, in what follows, together with the corresponding 
notation for rhombic patterns. 

Both packing problems find numerous counterparts in 
nature and in the laboratory; examples include biological 
microstructure 0j, the morphology of dry j5|| and wet 
foams in tubes, colloids in microchannels and the 
packing of fullerenes in nanotubes [1]. 

Simulations intended to search for the densest cylindri- 
cal arrangement of spheres Q have revealed a sequence 
of helical columnar phases, as the ratio R ^ D/d is var- 
ied (where d and D are the diameters of the hard spheres 
and the bounding cylinder, respectively). Some of these 
structures have been recognized in various experimental 
contexts, but there has until now been no clear under- 
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FIG. 1. (Color online), (a) Symmetric packing [5,4,1] - the 
black and red arrows indicate the periodicity V and primi- 
tive lattice vectors, respectively; (b) square packing [5,4]; (c) 
symmetric packing [9, 5, 4] which is connected to the structure 
[5, 4, 1] by an affine shear; (d) a line-slip structure between the 
symmetric packings [6,4,2] and [5,4,1]. 



standing of what determines the sequence of structural 
transitions as the cylinder diameter is increased. In this 
letter we demonstrate that the key to this progression is 
the familiar phyllotactic description. Using a straightfor- 
ward and elementary method we provide a full taxonomy 
of competing phases and also elucidate the geometric re- 
lationships between them. Surprisingly, we can use this 
method to also provide a full semiquantitative descrip- 
tion of the helical structures together with their packing 
densities. In general the densest structures are of the 
kind called "staggered" by Pickett [l]: we shall use the 
term "line slip". 
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We begin by discussing the disk packing problem. It 
is known that the densest packing arrangement of disks 
on a plane is that of a symmetric triangular lattice with 
lattice spacing d' (where d! is the disk diameter) [l^. A 
piece of the lattice can be excised and wrapped onto a 
cylinder of diameter D' (without any mismatch at the 
wrap edges) if there exist two lattice points separated by 
a vector V of magnitude, 

IVI^ttL*'. (1) 

Thus V plays the role of the periodicity vector, as shown 
in Fig ([1^) (where the lattice vertices correspond to the 
disk centers). The excised section is bounded by two par- 
allel lines perpendicular to V, these are the wrap edges 
and correspond to lines with 6 — Q and 6 — 27r, respec- 
tively, on the surface of a cylinder described by polar 

coordinates = ^,d,zj. 

In the phyllotactic scheme the periodicity vector V can 
be defined by a set of three, ordered, positive integers [l — 
m -\- n,m,rL), known as parastichy numbers in biology. 
Any pair of these numbers are the absolute values of the 
coefficients expressing V as a combination of two basis 
vectors, chosen from among three possibilities. Fig ([T^) 
shows an example for the structure [/,m, n] = [5,4,1] 
along with the basis vectors of the triangular lattice. 

The resulting structure, when wrapped onto a cylin- 
der, is a dense, homogenous packing - in the sense that 
all disk sites are equivalent - which we call a symmetric 
lattice packing. Between the discrete values of D', and 
fixed d' , for which such symmetric packings exist we first 
consider the most obvious adjustment to fit the cylin- 
der. In this the required periodicity is accommodated by 
an afHne transformation, or strain, of the lattice (more 
complex arrangements are possible but we do not find it 
necessary to consider them). The affine transformation 
distorts the equilateral triangle into an isosceles trian- 
gle by keeping the length of two adjacent sides fixed and 
varying the third. The imposed strain has three free pa- 
rameters; these can be chosen to satisfy Eq. {TJ and give 
each disk four contacting neighbors (in any one of three 
ways). 

As D' is varied, the strained structure proceeds from 
one symmetric packing to another, as shown in Fig ([T];). 
The rules for this process, when applied to the second 
and third phyllotactic indices, are as follows [ll| : 

(m, n) — !■ (m — n, n) 
(m, n) — > (m + n, m) 
(m, n) (m + n, n), 

where the new indices may have to be rearranged into de- 
scending order. For example, in the case of the symmet- 
ric packing [5,3,2] the above rules yield (3,2) — s> (2,1), 
(3,2) — > (5,3) and (3,2) — >• (5,2), where a rearrange- 
ment of the indices in the first case was necessary. Thus 



the symmetric packing [5,3,2] is connected, by an afBne 
shear, to the symmetric packings [3,2,1], [8,5,3], and 
[7,5,2]. 

The surface density (or area fraction of the disks on 
the plane) has a maximum value for a symmetric struc- 
ture. Between any two symmetric structures, connected 
in the way described above, the surface density has a 
lower value since each disk now has only four contacting 
neighbors. Such intermediate arrangements can be la- 
beled by using rhombic notation [p,q\ (see Q) and have 
a minimum surface density when the disks are arranged 
into a square lattice. Figure ^) shows such a square 
packing intermediate between two symmetric packings. 

We shall call these asymmetric lattice packings. In 
general, they are superseded by another type, as follows 
(indeed, one can show by convexity arguments that the 
asymmetric lattice structure is always unstable with re- 
spect to perturbations, under pressure). 

The second type of packing is achieved by an inhomo- 
geneous deformation of the symmetric lattice in which 
there is a localized strain or slip along a line (and its pe- 
riodic replicas) as shown in Fig ([TJi) . Again we have three 
choices of the direction of slip. The rest of the structure 
remains symmetric and close-packed, i.e. with six con- 
tacts for every disk. This line-slip distortion, as we shall 
call it, can again be expressed analytically. Other pos- 
sibilities, involving less localized strain exist, but these 
yield a lower density. Hence we include only the (opti- 
mized) case of single line-slip solutions. 

Again there is a simple rule for the close-packed struc- 
tures that are the endpoints of line-slip solutions; in view 
of their importance, we will identify explicitly the three 
types of solution |ll| : 

(1) (m, n) — )■ (m + 1, ri) or (to — l,r7,) 

(2) (m, n) — > (m, n + 1) or {m,n — l) 

(3) (m, n) — (m + 1, n — 1) or (m — l,n+l), 

where the leading numbers denote the direction ai, a.2 
or a3 of the line-slip. Again the above rules apply to the 
second and third phyllotactic indices of a given close- 
packed structure and keep either n, to or Z constant. For 
example, by using the above rules, the symmetric packing 
[5,3,2] is connected by a line slip along ai to [6,4,2] 
or [4,2,2], a line slip along a2 yields [4,3,1] or [6,3,3] 
and a line slip along as yields [5,4,1] or [5,3,2] (note 
in the second case along Sa that a rearrangement of the 
phyllotactic indices into descending order was necessary, 
the new structure has the same indices as the initial one) . 
Degeneracy is also possible; e.g. [4, 2, 2] is connected to 
[5, 3, 2] by a line slip along ai or a2. 

Such solutions (inspired by our simulations - see below) 
are the densest we have found by analytical searches, but 
we can offer no proof that excludes other denser struc- 
tures of greater complexity. Fig presents the ana- 
lytic results (to be published in detail elsewhere [ll| ) 
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FIG. 2. (Color online). Area fraction of disk packings on the 
plane which are consistent with wrapping onto a cylinder of 
diameter D' — \V\/-k. Asymmetric packings are shown by 
(red) dashed curves; line-slip solutions are (black) continuous 
curves. 



for the surface density of all of the packings described 
above, where we have set d! — 1. Asymmetric packings 
are shown by (red) dashed curves, line-slip solutions by 
(black) continuous curves; the symmetric packings are 
the labeled peaks with maximum density, while the lines 
connecting them are labeled by using the rhombic nota- 
tion. 

These results already show a strong resemblance to 
the corresponding simulation data for sphere packings 
in cylinders (but before we make a final comparison we 
will make a further transformation). The heavy contin- 
uous (black) line in Fig ([3]) indicates the computed op- 
timal packing density (or volume fraction) of monodis- 
perse hard spheres in an infinitely long cylindrical tube 
for R{= D/d). The present analysis was stimulated by 
that the numerical work of Pickett et al [1] , who used the 
method of simulated annealing to investigate the cylin- 
drical packing of hard spheres. We have employed two al- 
ternative approaches to verify and extend their findings. 
One is based on simulated annealing and employs the 
use of spiral boundary conditions to make the search for 
optimal structures more efficient, the other on sequential 
deposition. Both will be described elsewhere. The data 
of Fig ([3]) is obtained by the method with spiral bound- 
ary conditions, but the results of all three approaches are 
consistent. 

A representative sample of the simulated 3D structures 
is shown in Fig Below each structure is the "rolled 
out" phyllotactic pattern. The sphere centers lie on an in- 
ner cylinder D' = D — d and can be mapped to the plane 
by using Cartesian coordinates [{D'9)/2,z\, where they 
are shown as black dots. Contacting spheres are indi- 
cated by dots connected by black lines. The periodicity 




FIG. 3. (Color online). Comparison of simulated and an- 
alytic volume fractions for the densest packings, upper and 
lower curves respectively. Dotted lines are a guide for the 
eye to the detailed correspondence. In the analytic results, 
the line-slip structures are identified by black dashed lines 
and the red dashed-dotted lines denote asymmetric packing 
structures. The numbers 1,2 and 3 (see text) denote the type 
of line slip observed (upper curve) or predicted (lower curve) 
- with 1\2 denoting a degenerate case. 



vector, between two identical points on the cylinder, is 
denoted by a red arrow. The relationship to planar line- 
slip structures is evident. 

We relate the disk packings to those of spheres in con- 
tact with a cylinder, such as those shown in Fig 
as follows. Recall that the sphere centers lie on an in- 
ner cylinder of diameter D' . The sphere packing den- 
sity is trivially related to the point density of centers 
on that cylinder and, hence, to that of the planar pat- 
tern of points obtained by rolling out the cylinder on a 
plane. The qualitative resemblance to the disk problem 
becomes obvious, but its exact quantitative expression is 
not straightforward. Two sphere centers, separated by a 
vector r in 3D correspond to points in 2D whose separa- 
tion S depends on the orientation of r. For D' > d, the 
two extreme cases are S' = |r|, for r parallel to the cylin- 
der axis, and S — £)'sin~^(|r|/I?'), for r perpendicular 
to the cylinder axis. We adopt, as an ansatz, an inter- 
polation between these two limits which has the correct 
analytic form, 

[|r|sin(0)]V [D'sin-\\v\/D')coB{<t>)\\ (2) 

where (j) is the angle that the points on the inner cylin- 
der make with the vector V. This rule is applicable for 
R> 2 and a 2D packing that is consistent with it is sim- 
ply a packing of identical ellipses, whose major axes are 
all parallel to the vector V defined above, and may be 
obtained from any one of the disk packings by a simple 
stretch in the direction of V (which also changes V it- 
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FIG. 4. (Color online). Simulated 3D structures for various 
values of R, with the corresponding "rolled out" phyllotactic 
diagram below (the red arrow denotes the periodicity vector). 



self). In this way we arrive at a transformation of the 
circular disk packings into an approximate representa- 
tion of corresponding sphere packings, with densities as 
shown in Fig ([3]). 

This is a good approximation for our purposes; we note 
that the ratio between the predicted volume fraction, us- 
ing our method, and numerical results is already 0.9426 
at i? = 2.15 (i.e. the symmetric [3,3,0] structure) and 
continues to rapidly improve as |V| becomes larger. A 
more cumbersome method that uses the exact depen- 
dence of d on (p is perfectly possible, but a small sacrifice 
of accuracy preserves the extreme simplicity and trans- 
parency of our approach. By reference to the correspond- 
ing disk packing the approximate description developed 
above identifies the nature of all of the simulated sphere 
packing structures involved and is mostly correct in pick- 
ing the optimal structure. Transition between different 
sphere packings structures, indicated by the vertical dot- 
ted lines at the top of Fig can be understood by 
reference to the analogous disk packing problem. 

There are two situations leading to qualitative discrep- 
ancies between the theory and numerical results. The 
first corresponds to the region i? < 2.0 for which the ap- 
proximation Eq. 12]) fails. The result is that the twisted 
zigzag structure [9j arises unexpectedly as a stable asym- 
metric lattice structure. The transition to the line-slip 
structure at R = 1.995 is therefore of a different charac- 
ter from the others. Above this point, all of the inter- 
mediate structures are derived from the planar line-slip 
configurations. They are of the staggered spiral charac- 
ter (see the second and third figures in Fig dU), noted 
by Pickett et al in their simulations and exist up 
to R ~ 1 + 1/ sin(7r/5). Beyond this point not all the 
spheres are in contact with the boundary and the corre- 



spondence with the phyllotactic scheme is lost, although 
it may prove useful in the analysis of the surface layers. 

The second discrepancy is due to neglecting higher or- 
der corrections to S in Eq. ([2]) and leads to the obser- 
vation of line-slip solutions of type (2), in the numerical 
data, where type (3) would have been expected. This 
effect, whereby type (2) solutions have a greater than ex- 
pected volume fraction, diminishes for larger |V| as the 
higher order corrections to S fade away. 

A notable detail that emerges in the planar lattice 
analysis, and is retained in transformation to the sphere 
packing problem, is the singular behavior just before the 
points labeled 220, 330, etc. For reasons of symmetry, 
the decrease of density takes a square-root form and this 
is observed in the simulation data (most clearly in plots 
of the derivative of the density with respect to D). 

In summary, for2<i?<l-|-l/ sin(7r/5), each phyl- 
lotactic state (l,m,n) corresponds to a local density max- 
imum. Because of periodic constraints, and the nonover- 
lapping nature of hard spheres, these maxima are at- 
tained only at discrete values of D where six contacts are 
possible. Where only four contacts are possible, the line- 
slip structure with the highest density is observed. The 
density of the line-slip structure varies continuously with 
D, but transitions between them involve a structural dis- 
continuity. Though we are not aware of any experimental 
observations of line-slip structures, transitions between 
phyllotactic states have been observed in other systems 
|12| . In further work we will provide full details of the 
analytical results used here |ll| and include questions of 
stability and chirality. We shall extend our numerical re- 
sults to larger cylinder diameters - the packing problem 
then changes to one in which not all spheres touch the 
cylinder. 
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